** Graph comparing DPI and GWF

use DPI2012, clear
keep  countryname ifs year yrsoffc finittrm yrcurnt
gen cowcode = .
rename country country
do cowcodes
sum cow year
list country year if cow==.
drop if cow==.    /* South Sudan dropped */
drop if cow==352 /* Cyprus and Turk Cyprus */
recode cow (316=315) if year<1993   /* Czechoslakia code prior to 1993 */
recode cow (260=255)   /*   West German code */
tsset cow year
sort cow year
gen fail = F.yrsoffc==1 if yrsoffc~=-999
tab fail
gen DPI_fail = fail

merge cow year using GWFglobal
drop if year>2010 | year<1975    /* 1975-2010 are the years when the 2 data sets overlap */
drop if yrso==-999
tab _merge
browse if _merge==1   /* small countries not coded in GWF */
 
*Autocracy transition*
gen gwf_AA = gwf_fail==1 & gwf_next~="democracy" & gwf_next~="provisional"  & gwf_non~="democracy" & gwf_non~="provisional"
*Transition to Democracy*
gen gwf_AD = gwf_fail==1 & (gwf_next=="democracy" | gwf_next=="provisional")
*Democratic Failure*
gen gwf_DA = gwf_fail==1 & (gwf_next~="democracy" & gwf_next~="provisional")  & (gwf_non=="democracy" | gwf_non=="provisional")

tab gwf_fail gwf_AA
tab gwf_fail gwf_AD
tab gwf_fail gwf_DA


** Fix (4) inauguration date differences **
		*Panama 1990; Chile 1990; Uruguay 1985; South Korea 1988 
gen gwffail = gwf_fail
recode gwffail (0=1) if ((cowc==95 & year==1990) | (cowc==155 & year==1990) | (cowc==165 & year==1985) | (cowc==732 & year==1988)  )
recode gwf_AD (0=1) if  ((cowc==95 & year==1990) | (cowc==155 & year==1990) | (cowc==165 & year==1985) | (cowc==732 & year==1988)  )
recode gwf_AD (1=0) if  ((cowc==95 & year==1989) | (cowc==155 & year==1989) | (cowc==165 & year==1984) | (cowc==732 & year==1987)  )

**leader failure years when NO regime collapse**
tab DPI_fail if gwffail==0 & gwf_regimetype~="NA"   & DPI_fail~=.
tab DPI_fail if gwffail==0 & (gwf_non=="democracy" | gwf_non=="provisional")  & DPI_fail~=.

**leader failure years when YES regime collapse**
tab DPI_fail if gwffail==1 & gwf_AA==1 & gwf_regimetype~="NA" & DPI_fail~=.
tab DPI_fail if gwffail==1 & gwf_AD==1 & gwf_regimetype~="NA" & DPI_fail~=.
tab DPI_fail if gwffail==1 & gwf_DA==1 & (gwf_non=="democracy" | gwf_non=="provisional") & DPI_fail~=.

gen graphgwfdpi= .
replace graphgwfdpi = 428 in 1    /* Democracy survives */
replace graphgwfdpi  = . in 2    
replace graphgwfdpi  = 22 in 3     /* Democracy-Autocracy */
replace graphgwfdpi = . in 4    
replace graphgwfdpi  = 55 in 5     /* Autocracy-Democracy */
replace graphgwfdpi  = . in 6  
replace graphgwfdpi  = 58 in 7    /* Autocracy-Autocracy */
replace graphgwfdpi  = . in 8  
replace graphgwfdpi  = 171 in 9    /* Autocracy survives */

gen index  = _n in 1/9
label define allfails  1 "Democracy survives" 2 " "  3 "Democracy-Autocracy"  4 " " 5 "Autocracy-Democracy" /*
*/ 6 " " 7 "Autocracy-Autocracy" 8 " "  9 "Autocracy survives"   
label values index  allfails 

twoway (bar graphgwfdpi  index, scheme(lean2)  ylabel(0(100)400,glcolor(gs14))   /*
*/ xlabel(1(1)9, valuelabel labsize(small) labcolor(black) labgap(tiny) noticks) xsize(10) ysize(5) )/*
*/ (scatter graphgwfdpi  index, ms(none) mla(graphgwfdpi ) mlabpos(12) graphregion(color(white)) /*
*/ xtitle("Type of event",height(8)) ytitle("Number of country-years")  xscale(range (0 10)) yscale(range(0 450)) /*
*/ legend(pos(12) col(1) ring(0) label(1 "DPI Government Failure") label(2 "") bmargin(0)) )
*graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\DPIGWF.pdf", as(pdf)  replace




*****************************************************
** Merge and clean data for Ahmed 2012 replication **
*****************************************************
use AhmedAPSR, clear
set more off
rename COUNTRY_NAME country
qui gen cowcode=.
qui do cowcodes
sum year cow
qui egen tag = tag(country) if cow==.
list country if tag==1 , clean
qui drop if cow==.
sum cow year
sort cow year
merge cow year using GWFglobal
qui dprobit turnover aid_remit beta beta_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy*, cluster(govtcode)
gen sumsamp=e(sample)
tab _merge if sumsamp  /*Need to get Ethiopia post 1993*/

tab turnover gwf_fail if sumsamp==1 & turnover==1, m
/*turnover with disagreedpi==. are small island countries*/
gen turnover1 = turnover
replace turnover1 = 0 if gwf_fail==1   
gen turnover2 = turnover
replace turnover2 = 0 if gwf_fail==0 | (gwf_fail==. & turnover==1) 


/*One-year election/national conference differences*/
recode turnover1 (1=0) if turnover==1 & ((cow==339 & year==1992) | (cow==434 & year==1991 ) | (cow==484 & year==1992) | (cow==436 & year==1993 ))
recode turnover2 (0=1) if turnover==1 & ((cow==339 & year==1992) | (cow==434 & year==1991 ) | (cow==484 & year==1992) | (cow==436 & year==1993 ))

/*Cote d'Ivoire error: Ahmed codes Cote d'Ivoire 1984 as turnover but this isn't correct and isn't in the original DPI data */
replace turnover2 = 0 if cow==437 & year==1984  
replace turnover1 = 1 if cow==437 & year==1984

tsset cow year
qui dprobit turnover aid_remit beta beta_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy*, cluster(govtcode)
gen lpolity = l.POLITY2
gen beta1 = (((POLITY2+11)*-1) +22)/21   /*DEC 31 of obs year coding of polity*/
gen beta2 = (((lpolity+11)*-1) +22)/21   /*JAN 1 of obs year coding of polity*/
gen beta3 = l.beta                       /*JAN 1 of obs year coding of polity, for transformed index */
sum beta POLITY2
label var beta "Autocracy score"
label var POLITY2 "Polity index"
twoway scatter beta POLITY2
*graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\AutocracyScore.pdf", as(pdf)                            replace

gen beta1_AR = beta1*aid_remit
gen beta2_AR = beta2*aid_remit
gen beta3_AR = beta3*aid_remit

corr beta beta1 beta2 beta3
corr beta_AR beta1_AR beta2_AR beta3_AR

gen beta_muslimoil = beta*muslimoil
gen beta1_muslimoil = beta1*muslimoil
gen beta2_muslimoil = beta2*muslimoil
gen beta3_muslimoil = beta3*muslimoil

label var aid_remit "Aid & Remittances"
label var beta3 "Autocracy"
label var turnover "Failure"
label var turnover1 "Failure"
label var turnover2 "Failure"
drop _merge
saveold temp, replace
saveold temp_Ahmed, replace


**********************************
** Ahmed 2012, Table 3 column 2 **
**********************************

sort country year
browse country year if sumsamp==1 & turnover1==1 & turnover==1
browse country year if sumsamp==1 & turnover2==1 & turnover==1

*Replication
   qui probit turnover aid_remit beta beta_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy* if sumsamp==1, cluster(govtcode)
	*sum beta if e(sample), detail
	tab turnover if e(sample)
	lincom aid_remit + beta_AR*.048   /*5 pctile*/
	lincom aid_remit + beta_AR*.5   /*95 pctile*/
	est store a1

*DPI fail and GWF fail
  qui probit turnover2 aid_remit beta beta_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy* if sumsamp==1, cluster(govtcode)
	*sum beta if e(sample), detail
	tab turnover2 if e(sample)
	lincom aid_remit + beta_AR*.053   /*5 pctile*/
	lincom aid_remit + beta_AR*.5   /*95 pctile*/
	est store a2

*DPI fail but NO GWF fail
  qui probit turnover1 aid_remit beta beta_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy* if sumsamp==1, cluster(govtcode)
	*sum beta if e(sample), detail
	tab turnover1 if e(sample)
	lincom aid_remit + beta_AR*.048   /*5 pctile*/
	lincom aid_remit + beta_AR*.5   /*95 pctile*/
	est store a3

*Replication w. correctly lagged beta
  qui probit turnover aid_remit beta3 beta3_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy* if sumsamp==1, cluster(govtcode)
	*sum beta3 if e(sample), detail
	tab turnover if e(sample)
	lincom aid_remit + beta3_AR*.048   /*5 pctile*/
	lincom aid_remit + beta3_AR*.5   /*95 pctile*/
	est store a4

*DPI fail and GWF fail w. correctly lagged beta
  qui probit turnover2 aid_remit beta3 beta3_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy* if sumsamp==1, cluster(govtcode)
	*sum beta3 if e(sample), detail
	tab turnover2 if e(sample)
	lincom aid_remit + beta3_AR*.053   /*5 pctile*/
	lincom aid_remit + beta3_AR*.5   /*95 pctile*/
	est store a5
 
*DPI fail and NO GWF fail w. correctly lagged beta
  qui probit turnover1 aid_remit beta3 beta3_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy* if sumsamp==1, cluster(govtcode)
	*sum beta3 if e(sample), detail
	tab turnover1 if e(sample)
	lincom aid_remit + beta3_AR*.048   /*5 pctile*/
	lincom aid_remit + beta3_AR*.5   /*95 pctile*/
	est store a6
 
set matsize 2000
*set emptycells drop
estout a1 a2 a3 a4 a5 a6 using Table3.tex, cells(b(star  fmt(%9.3f)) se(par fmt(%9.2f))) stats(ll r2 N) style(tex) replace label starlevels(* 0.10 ** 0.05)
set matsize 500

********************************
** Full sample specifications **
********************************
use temp, clear
qui probit turnover aid_remit beta3 beta3_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time time2 time3 year  africa nam asia sam if sumsamp==1, cluster(cow)
sum beta3 if e(sample), detail
lincom aid_remit + beta3_AR*.048  /*5 pctile*/
lincom aid_remit + beta3_AR*.5   /*95 pctile*/
 
qui probit turnover2 aid_remit beta3 beta3_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time time2 time3 year  africa nam asia sam if sumsamp==1, cluster(cow)
sum beta3 if e(sample), detail
lincom aid_remit + beta3_AR*.053   /*5 pctile*/
lincom aid_remit + beta3_AR*.5   /*95 pctile*/
 
qui probit turnover1 aid_remit beta3 beta3_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time time2 time3 year  africa nam asia sam if sumsamp==1, cluster(cow)
sum beta3 if e(sample), detail
lincom aid_remit + beta3_AR*.048   /*5 pctile*/
lincom aid_remit + beta3_AR*.5   /*95 pctile*/
 
 *******************
 *** Simulations ***   ydummy53, codummy28 for turnover2, codummy54 for turnover/turnover1 time_dum6
 *******************
 /*    Simulations run separately in Stata 9  
 
 * (1) *
 use temp, clear
 set more off
 probit turnover aid_remit beta beta_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy* if sumsamp==1, cluster(govtcode)
 matrix m1 = e(b)
 scalar s1 = colsof(m1)
 scalar list s1 /*162*/
 set seed 9879789
 drawnorm b1-b162, n(10000) means(e(b)) cov(e(V))  clear 
 gen PROB1dem=.
 gen PROB1demlow=.
 gen PROB1demhigh=.
 gen PROB1aut=.
 gen PROB1autlow=.
 gen PROB1authigh=.

 gen a =.
 gen ar_axis = (_n-1)/5 + 0.5 if _n<=132
 /*We want nt to run from 0.5 to 26.5, ie from 5pctile to 95pctile*/ 
 local a=0.5
 while `a' <= 26.5 {
  	gen x_betaDem  = b1*`a'+b2*.05+b3*.05*`a'+b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b15*1+b61*1+b156*1+b162*1
  	gen x_betaAut  = b1*`a'+b2*.5+b3*.5*`a'  +b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b15*1+b61*1+b156*1+b162*1
 	gen probD = normal(x_betaDem)
 	gen probA = normal(x_betaAut) 
 	egen probDhat = mean(probD)
      _pctile probD, p(2.5,97.5) 
      scalar probDlow= r(r1) 
      scalar probDhigh= r(r2) 
 	egen probAhat = mean(probA)
      _pctile probA, p(2.5,97.5) 
      scalar probAlow= r(r1) 
      scalar probAhigh= r(r2) 
   	quietly replace a = `a' 
 	quietly replace PROB1dem =  probDhat  if ar_axis==a
 	quietly replace PROB1demlow =  probDlow  if ar_axis==a
 	quietly replace PROB1demhigh =  probDhigh  if ar_axis==a
 	quietly replace PROB1aut =  probAhat  if ar_axis==a
 	quietly replace PROB1autlow =  probAlow  if ar_axis==a
 	quietly replace PROB1authigh =  probAhigh  if ar_axis==a
 	drop x_beta* prob*  
 	local a = `a' + .2
 }
 label var ar_axis "Aid + Remit"
 label var PROB1dem "All failures"
 label var PROB1aut "All failures"
 sort ar_axis
 twoway (line PROB1dem ar_axis)  (line PROB1aut ar_axis) 
 saveold Arep1, replace
 
  * (2) *
  use temp, clear
  set more off
  probit turnover2 aid_remit beta beta_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis   time_dum*  codummy* ydummy* if sumsamp==1, cluster(govtcode)
   matrix m1 = e(b)
  scalar s1 = colsof(m1)
  scalar list s1 /*110*/
  drawnorm b1-b110, n(10000) means(e(b)) cov(e(V))  clear 
 gen PROB2dem=.
 gen PROB2demlow=.
 gen PROB2demhigh=.
 gen PROB2aut=.
 gen PROB2autlow=.
 gen PROB2authigh=.
  gen a =.
  gen ar_axis = (_n-1)/5 + 0.5 if _n<=132
  /*We want nt to run from 0.5 to 26.5, ie from 5pctile to 95pctile*/ 
  local a=0.5
  while `a' <= 26.5 {
  	gen x_betaDem  = b1*`a'+b2*.05+b3*.05*`a'+b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b15*1+b42*1+b104*1+b110*1
  	gen x_betaAut  = b1*`a'+b2*.5+b3*.5*`a'  +b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b15*1+b42*1+b104*1+b110*1
  	gen probD = normal(x_betaDem)
  	gen probA = normal(x_betaAut) 
 	egen probDhat = mean(probD)
      _pctile probD, p(2.5,97.5) 
      scalar probDlow= r(r1) 
      scalar probDhigh= r(r2) 
 	egen probAhat = mean(probA)
      _pctile probA, p(2.5,97.5) 
      scalar probAlow= r(r1) 
      scalar probAhigh= r(r2) 
   	quietly replace a = `a' 
 	quietly replace PROB2dem =  probDhat  if ar_axis==a
 	quietly replace PROB2demlow =  probDlow  if ar_axis==a
 	quietly replace PROB2demhigh =  probDhigh  if ar_axis==a
 	quietly replace PROB2aut =  probAhat  if ar_axis==a
 	quietly replace PROB2autlow =  probAlow  if ar_axis==a
 	quietly replace PROB2authigh =  probAhigh  if ar_axis==a
  	drop x_beta* prob*  
  	local a = `a' + .2
  }
  label var ar_axis "Aid + Remit"
  label var PROB2dem "Regime collapse"
  label var PROB2aut "Regime collapse"
  sort ar_axis
  twoway (line PROB2dem ar_axis)  (line PROB2aut ar_axis) 
 saveold Arep2, replace
 
  * (3) *
  use temp, clear
  set more off
  probit turnover1 aid_remit beta beta_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum*  codummy* ydummy* if sumsamp==1, cluster(govtcode)
    matrix m1 = e(b)
  scalar s1 = colsof(m1)
  scalar list s1 /*145*/
  drawnorm b1-b145, n(10000) means(e(b)) cov(e(V))  clear 
  gen PROB3dem=.
  gen PROB3demlow=.
  gen PROB3demhigh=.
  gen PROB3aut=.
  gen PROB3autlow=.
  gen PROB3authigh=.
  gen a =.
  gen ar_axis = (_n-1)/5 + 0.5 if _n<=132
  /*We want nt to run from 0.5 to 26.5, ie from 5pctile to 95pctile*/ 
  local a=0.5
  while `a' <= 26.5 {
  	gen x_betaDem  = b1*`a'+b2*.05+b3*.05*`a'+b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b15*1+b54*1+b140*1+b145*1
  	gen x_betaAut  = b1*`a'+b2*.5+b3*.5*`a'  +b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b15*1+b54*1+b140*1+b145*1
  	gen probD = normal(x_betaDem)
  	gen probA = normal(x_betaAut) 
 	egen probDhat = mean(probD)
      _pctile probD, p(2.5,97.5) 
      scalar probDlow= r(r1) 
      scalar probDhigh= r(r2) 
 	egen probAhat = mean(probA)
      _pctile probA, p(2.5,97.5) 
      scalar probAlow= r(r1) 
      scalar probAhigh= r(r2) 
   	quietly replace a = `a' 
 	quietly replace PROB3dem =  probDhat  if ar_axis==a
 	quietly replace PROB3demlow =  probDlow  if ar_axis==a
 	quietly replace PROB3demhigh =  probDhigh  if ar_axis==a
 	quietly replace PROB3aut =  probAhat  if ar_axis==a
 	quietly replace PROB3autlow =  probAlow  if ar_axis==a
 	quietly replace PROB3authigh =  probAhigh  if ar_axis==a
  	drop x_beta* prob*  
  	local a = `a' + .2
  }
  label var ar_axis "Aid + Remit"
  label var PROB3dem "No regime collapse"
  label var PROB3aut "No regime collapse"
  sort ar_axis
  twoway (line PROB3dem ar_axis)  (line PROB3aut ar_axis) 
  saveold Arep3, replace
 
 
  * (4) *
 use temp, clear
 set more off
 probit turnover aid_remit beta3 beta3_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum*  codummy* ydummy*  if sumsamp==1, cluster(govtcode)
 drawnorm b1-b162, n(10000) means(e(b)) cov(e(V))  clear 
 gen PROB4dem=.
 gen PROB4demlow=.
 gen PROB4demhigh=.
 gen PROB4aut=.
 gen PROB4autlow=.
 gen PROB4authigh=.
 gen a =.
 gen ar_axis = (_n-1)/5 + 0.5 if _n<=132
 /*We want nt to run from 0.5 to 26.5, ie from 5pctile to 95pctile*/ 
 local a=0.5
 while `a' <= 26.5 {
  	gen x_betaDem  = b1*`a'+b2*.05+b3*.05*`a'+b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b15*1+b61*1+b156*1+b162*1
  	gen x_betaAut  = b1*`a'+b2*.5+b3*.5*`a'  +b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b15*1+b61*1+b156*1+b162*1
 	gen probD = normal(x_betaDem)
 	gen probA = normal(x_betaAut) 
 	egen probDhat = mean(probD)
      _pctile probD, p(2.5,97.5) 
      scalar probDlow= r(r1) 
      scalar probDhigh= r(r2) 
 	egen probAhat = mean(probA)
      _pctile probA, p(2.5,97.5) 
      scalar probAlow= r(r1) 
      scalar probAhigh= r(r2) 
   	quietly replace a = `a' 
 	quietly replace PROB4dem =  probDhat  if ar_axis==a
 	quietly replace PROB4demlow =  probDlow  if ar_axis==a
 	quietly replace PROB4demhigh =  probDhigh  if ar_axis==a
 	quietly replace PROB4aut =  probAhat  if ar_axis==a
 	quietly replace PROB4autlow =  probAlow  if ar_axis==a
 	quietly replace PROB4authigh =  probAhigh  if ar_axis==a
 	drop x_beta* prob*  
 	local a = `a' + .2
 }
 label var ar_axis "Aid + Remit"
 label var PROB4dem "All failures"
 label var PROB4aut "All failures"
 sort ar_axis
 twoway (line PROB4dem ar_axis)  (line PROB4aut ar_axis) 
 saveold Arep4, replace
 
 * (5) *
 use temp, clear
 set more off
 probit turnover2 aid_remit beta3 beta3_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum*  codummy* ydummy* if sumsamp==1, cluster(govtcode)
   matrix m1 = e(b)
 scalar s1 = colsof(m1)
 scalar list s1 /*109*/
 drawnorm b1-b109, n(10000) means(e(b)) cov(e(V))  clear 
 gen PROB5dem=.
 gen PROB5demlow=.
 gen PROB5demhigh=.
 gen PROB5aut=.
 gen PROB5autlow=.
 gen PROB5authigh=.
 gen a =.
 gen ar_axis = (_n-1)/5 + 0.5 if _n<=132
 /*We want nt to run from 0.5 to 26.5, ie from 5pctile to 95pctile*/ 
 local a=0.5
 while `a' <= 26.5 {
  	gen x_betaDem  = b1*`a'+b2*.05+b3*.05*`a'+b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b15*1+b41*1+b105*1+b109*1
  	gen x_betaAut  = b1*`a'+b2*.5+b3*.5*`a'  +b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b15*1+b41*1+b105*1+b109*1
 	gen probD = normal(x_betaDem)
 	gen probA = normal(x_betaAut) 
 	egen probDhat = mean(probD)
      _pctile probD, p(2.5,97.5) 
      scalar probDlow= r(r1) 
      scalar probDhigh= r(r2) 
 	egen probAhat = mean(probA)
      _pctile probA, p(2.5,97.5) 
      scalar probAlow= r(r1) 
      scalar probAhigh= r(r2) 
   	quietly replace a = `a' 
 	quietly replace PROB5dem =  probDhat  if ar_axis==a
 	quietly replace PROB5demlow =  probDlow  if ar_axis==a
 	quietly replace PROB5demhigh =  probDhigh  if ar_axis==a
 	quietly replace PROB5aut =  probAhat  if ar_axis==a
 	quietly replace PROB5autlow =  probAlow  if ar_axis==a
 	quietly replace PROB5authigh =  probAhigh  if ar_axis==a
 	drop x_beta* prob*  
 	local a = `a' + .2
 }
 label var ar_axis "Aid + Remit"
 label var PROB5dem "Regime collapse"
 label var PROB5aut "Regime collapse"
 sort ar_axis
 twoway (line PROB5dem ar_axis)  (line PROB5aut ar_axis) 
 saveold Arep5, replace
 
  * (6) *
 use temp, clear
 set more off
 probit turnover1 aid_remit beta3 beta3_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum*  codummy* ydummy*  if sumsamp==1, cluster(govtcode)
 drawnorm b1-b145, n(10000) means(e(b)) cov(e(V))  clear 
 gen PROB6dem=.
 gen PROB6demlow=.
 gen PROB6demhigh=.
 gen PROB6aut=.
 gen PROB6autlow=.
 gen PROB6authigh=.
 gen a =.
 gen ar_axis = (_n-1)/5 + 0.5 if _n<=132
 /*We want nt to run from 0.5 to 26.5, ie from 5pctile to 95pctile*/ 
 local a=0.5
 while `a' <= 26.5 {
  	gen x_betaDem  = b1*`a'+b2*.05+b3*.05*`a'+b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b15*1+b54*1+b140*1+b145*1
  	gen x_betaAut  = b1*`a'+b2*.5+b3*.5*`a'  +b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b15*1+b54*1+b140*1+b145*1
 	gen probD = normal(x_betaDem)
 	gen probA = normal(x_betaAut) 
 	egen probDhat = mean(probD)
      _pctile probD, p(2.5,97.5) 
      scalar probDlow= r(r1) 
      scalar probDhigh= r(r2) 
 	egen probAhat = mean(probA)
      _pctile probA, p(2.5,97.5) 
      scalar probAlow= r(r1) 
      scalar probAhigh= r(r2) 
   	quietly replace a = `a' 
 	quietly replace PROB6dem =  probDhat  if ar_axis==a
 	quietly replace PROB6demlow =  probDlow  if ar_axis==a
 	quietly replace PROB6demhigh =  probDhigh  if ar_axis==a
 	quietly replace PROB6aut =  probAhat  if ar_axis==a
 	quietly replace PROB6autlow =  probAlow  if ar_axis==a
 	quietly replace PROB6authigh =  probAhigh  if ar_axis==a
 	drop x_beta* prob*  
 	local a = `a' + .2
 }
 label var ar_axis "Aid + Remit"
 label var PROB6dem "No regime collapse"
 label var PROB6aut "No regime collapse"
 sort ar_axis
 twoway (line PROB6dem ar_axis)  (line PROB6aut ar_axis) 
 saveold Arep6, replace
 */
 
 use Arep1, clear
 sort ar_axis
 merge ar_axis using Arep2
 sort ar_axis
 drop _merge
 merge ar_axis using Arep3
 sort ar_axis
 drop _merge
 merge ar_axis using Arep4
 sort ar_axis
 drop _merge
 merge ar_axis using Arep5
 sort ar_axis
 drop _merge
 merge ar_axis using Arep6
 sort ar_axis
 drop _merge
 merge using temp
 set scheme lean1
 label var aid_remit "Aid + Remit distribution"
 sort ar_axis
 drop _merge
 
twoway  (hist aid_remit if aid_remit<=26.5 & sumsam==1 & aid_remit>=0.5, lcolor(gs12) bin(100)) /*
*/ (line PROB1aut ar_axis if ar_axis<=26.5, yaxis(1))  /*
*/ (line PROB2aut ar_axis if ar_axis<=26.5, yaxis(1))  /*
*/ (line PROB3aut ar_axis if ar_axis<=26.5, yaxis(1))  /*
*/ ,  xtitle("Aid + Remit (%GDP)", size(3.5)) ytitle("Pr(Failure)", size(3.5)) /*
*/  legend(pos(12) col(2) ring(1) label(1 "Aid + Remit distribution")  label(2 "All failures") label(3 "Regime collapse") label(4 "No regime collapse"))  xscale(range (0 26)) xlabel(0 (5) 25)
**graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\AhmedReplication.pdf", as(pdf)                            replace

twoway  (hist aid_remit if aid_remit<=26.5 & sumsam==1 & aid_remit>=0.5, lcolor(gs12) bin(100)) /*
*/ (line PROB4aut ar_axis if ar_axis<=26.5, yaxis(1))  /*
*/ (line PROB5aut ar_axis if ar_axis<=26.5, yaxis(1))  /*
*/ (line PROB6aut ar_axis if ar_axis<=26.5, yaxis(1))  /*
*/ ,  xtitle("Aid + Remit (%GDP)", size(3.5)) ytitle("Pr(Failure)", size(3.5)) /*
*/  legend(pos(12) col(2) ring(1) label(1 "Aid + Remit distribution")  label(2 "All failures") label(3 "Regime collapse") label(4 "No regime collapse"))  xscale(range (0 26)) xlabel(0 (5) 25)
**graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\AhmedCorrected.pdf", as(pdf)                            replace



******************************
** 95% Confidence Intervals **
******************************

twoway  (hist aid_remit if aid_remit<=26.5 & sumsam==1 & aid_remit>=0.5, lcolor(gs12) bin(100)) /*
*/ (line PROB1aut PROB1autlow PROB1authigh ar_axis if ar_axis<=26.5, clpattern(solid dash dash) yaxis(1)) /*
*/ ,  xtitle("Aid + Remit (%GDP)", size(4)) ytitle("Pr(Failure)", size(4)) /*
*/ legend(off)
*graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\AhmedReplication_all95CI.pdf", as(pdf)                            replace

twoway  (hist aid_remit if aid_remit<=26.5 & sumsam==1 & aid_remit>=0.5, lcolor(gs12) bin(100)) /*
*/ (line PROB2aut PROB2autlow PROB2authigh ar_axis if ar_axis<=26.5, clpattern(solid dash dash) yaxis(1)) /*
*/ ,  xtitle("Aid + Remit (%GDP)", size(4)) ytitle("Pr(Failure)", size(4)) /*
*/ legend(off)
*graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\AhmedReplication_fail95CI.pdf", as(pdf)                            replace

twoway  (hist aid_remit if aid_remit<=26.5 & sumsam==1 & aid_remit>=0.5, lcolor(gs12) bin(100)) /*
*/ (line PROB3aut PROB3autlow PROB3authigh ar_axis if ar_axis<=26.5, clpattern(solid dash dash) yaxis(1)) /*
*/ ,  xtitle("Aid + Remit (%GDP)", size(4)) ytitle("Pr(Failure)", size(4)) /*
*/ legend(off)
*graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\AhmedReplication_nofail95CI.pdf", as(pdf)                            replace

twoway  (hist aid_remit if aid_remit<=26.5 & sumsam==1 & aid_remit>=0.5, lcolor(gs12) bin(100)) /*
*/ (line PROB4aut PROB4autlow PROB4authigh ar_axis if ar_axis<=26.5, clpattern(solid dash dash) yaxis(1)) /*
*/ ,  xtitle("Aid + Remit (%GDP)", size(4)) ytitle("Pr(Failure)", size(4)) /*
*/ legend(off)
*graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\AhmedCorrected_all95CI.pdf", as(pdf)                            replace

twoway  (hist aid_remit if aid_remit<=26.5 & sumsam==1 & aid_remit>=0.5, lcolor(gs12) bin(100)) /*
*/ (line PROB5aut PROB5autlow PROB5authigh ar_axis if ar_axis<=26.5, clpattern(solid dash dash) yaxis(1)) /*
*/ ,  xtitle("Aid + Remit (%GDP)", size(4)) ytitle("Pr(Failure)", size(4)) /*
*/ legend(off)
*graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\AhmedCorrected_fail95CI.pdf", as(pdf)                            replace

twoway  (hist aid_remit if aid_remit<=26.5 & sumsam==1 & aid_remit>=0.5, lcolor(gs12) bin(100)) /*
*/ (line PROB6aut PROB6autlow PROB6authigh ar_axis if ar_axis<=26.5, clpattern(solid dash dash) yaxis(1)) /*
*/ ,  xtitle("Aid + Remit (%GDP)", size(4)) ytitle("Pr(Failure)", size(4)) /*
*/ legend(off)
*graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\AhmedCorrected_nofail95CI.pdf", as(pdf)                            replace


 
**********************************
** Ahmed 2012, Table 3 column 3 **
**********************************
 /* done separately in Stata 9 due to convergence issues; ivprobit are pretty unstable with two endogenous variables/ better to use linear...
use temp, clear
global c = "finittrm lngdpcap ggdpcap lnpop ldis hdis time time2 time3 ydummy* asia africa nam sam"
global p = " sam ggdpcap time time2 lnpop lngdpcap war africa "

 **Replication model
 qui ivprobit turnover (beta_AR=muslimoil) $c if sumsamp==1 , cluster(govtcode)
 lincom beta_AR
 est store c1
 
 **Use correctly lagged Autocracy score
 qui ivprobit turnover (beta3_AR=muslimoil) $c if sumsamp==1 , cluster(govtcode)
 lincom beta3_AR
 est store c2
 
 **Use correctly lagged Autocracy score  and parsimonious specification
 	* Drops: time3, finitterm, asia, nam, ydummy*, ldis, and hdis
 qui ivprobit turnover (beta3_AR   =muslimoil)  $p     if sumsamp==1 , cluster(govtcode) difficult technique(bfgs)
 lincom beta3_AR
 est store c3
 
 **Use correctly lagged Autocracy score  and parsimonious specification and add constituent terms
 	* Drops: time3, finitterm, asia, nam, ydummy*, ldis, and hdis 
 qui ivprobit turnover (beta3_AR aid_remit  =muslimoil beta3_muslim) beta3 $p   if sumsamp==1 , cluster(govtcode) difficult technique(bfgs)
 lincom beta3_AR
 lincom aid_remit
 lincom aid_remit + beta3_AR*0.048
 lincom aid_remit + beta3_AR*0.5
 est store c4
 
 **DPI fail and YES GWF fail: Use correctly lagged Autocracy score  and parsimonious specification and add constituent terms
 qui ivprobit turnover2 (beta3_AR aid_remit  =muslimoil beta3_muslim)   beta3  $p    if sumsamp==1 , cluster(govtcode) difficult technique(bfgs)
 lincom beta3_AR
 lincom aid_remit
 lincom aid_remit + beta3_AR*0.048
 lincom aid_remit + beta3_AR*0.5
 est store c5
 
 **DPI fail and NO GWF fail: Use correctly lagged Autocracy score  and parsimonious specification and add constituent terms
 qui ivprobit turnover1 (beta3_AR aid_remit  =muslimoil beta3_muslim) beta3 $p    if sumsamp==1 , cluster(govtcode) difficult technique(bfgs)
 lincom beta3_AR
 lincom aid_remit
 lincom aid_remit + beta3_AR*0.048
 lincom aid_remit + beta3_AR*0.5
 est store c6
estout c1 c2 c3 c4 c5 c6 using TableC1.tex, cells(b(star  fmt(%9.3f)) se(par fmt(%9.2f))) stats(ll r2 N) style(tex) replace label starlevels(* 0.10 ** 0.05)
*/

 *Double-check two-stage in a linear probability model: No results one way or another
 global c = "finittrm lngdpcap ggdpcap lnpop ldis hdis time time2 time3 ydummy* asia africa nam sam"

 qui ivreg2 turnover (aid_remit=muslimoil) $c if sumsamp==1, cluster(govtcode)    
 lincom aid_remit
 
 qui ivreg2 turnover (beta3_AR aid_remit=muslimoil beta3_muslimoil) beta3_AR $c if sumsamp==1, cluster(govtcode) 
 lincom aid_remit + beta3_AR*.5
 lincom aid_remit + beta3_AR*1
 
 qui ivreg2 turnover2 (beta3_AR aid_remit=muslimoil beta3_muslimoil) beta3_AR $c if sumsamp==1, cluster(govtcode)   
 lincom aid_remit + beta3_AR*.5
 lincom aid_remit + beta3_AR*1
 
 qui ivreg2 turnover1 (beta3_AR aid_remit=muslimoil beta3_muslimoil) beta3_AR $c if sumsamp==1, cluster(govtcode)  
 lincom aid_remit + beta3_AR*.5
lincom aid_remit + beta3_AR*1
